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Abstract 

We present an algorithm for measurement of the Green's function in the hy- 
bridization expansion continuous-time quantum Monte-Carlo based on con- 
tinuous estimators. Compared to the standard method, the present algo- 
rithm has similar or better accuracy with improvement notable especially 
at low perturbation orders and high frequencies. The resulting statistical 
noise is weakly correlated so high-accuracy data for the numerical analytical 
continuation can be produced. 

Keywords: CT-HYB, Measurement 



1. Introduction 

The Monte- Carlo simulations are among the most widely used techniques 
for investigation of quantum impurity models. The progress of the dynam- 
ical mean-field theory [l], [i], [sj in the past decade, which allowed material 
specific studies of multi-band Hubbard models, is to a large extent governed 
by the availability of reliable and fast solvers for multi-orbital impurity prob- 
lems. New generation quantum Monte-Carlo (QMC) techniques based on the 
stochastic sampling of perturbative expansions 0, lEi |6[ 0] quickly became the 
standard of quantum impurity simulations. 

Several QMC algorithms are currently used in this context differing by the 
type of expansion and symmetry of the impurity Hamiltonian. As the sim- 
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ulation of impurities with arbitrary interaction js], lof is computationally de- 
manding, specialized algorithms that can handle a simplified density-density 
interaction at much lower computational cost still play an important role. 

In this article we describe an improved technique for measurement of 
Green's functions in the strong-coupling CT-QMC algorithm (also called 
CT-HYB) for density-density interactions. The development was motivated 
by poor performance of the standard measurement technique in systems with 
strongly imbalanced perturbation orders for different orbitals. Such a situa- 
tion is common in materials containing orbitals with fluctuating occupation 
coupled by interaction to filled or empty orbitals, e.g. the tog orbitals in nick- 



elates [10|], Cg orbitals in early transition metal oxides [ll|, ll2|, t2g orbitals 



in the low-spin state of LaCoOa 
splitting ll 
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13|, or Hubbard model with crystal-field 



Although the filled orbitals experience little fluctuations 
themselves and their influence on the active (partially filled) orbitals is well 
approximated by a static potential, it is their Green's functions that may 
become prohibitively noisy. In some cases one can simply remove the filled 
orbitals from the effective model, e.g. early transition metals are commonly 
studied with t2g-only models [l6|. Often, however, some orbitals became 
filled/empty only in a certain part of the phase diagram and then they have 
to be kept in the model, e.g. to study spin/orbital ordering, or because 
one is interested in their dynamics, which may exhibit non-trivial dynamical 



renormalization 17, 10 



The paper is organized as follows. In section [2] we review the basic of the 
CT-HYB algorithm and the standard measurement technique. In section [3] 
we discuss on a general level the relationship between the expansions for the 
partition function and the Green's function. The new measurement algorithm 
is described in section H] with details and implementation notes summarized 
in section 13 Finally, in section E] we provide examples and performance tests. 



2. CT-HYB 

At this point, we briefly review the CT-HYB algorithm introduced in 
0, [sj (for details see comprehensive review 

We wish to study the dynamics of the (i-electrons in the Anderson impu- 
rity model with density- density interaction described by the Hamiltonian 

H ^ Y,(^aknak + Y,na(^E^ + Y,Ual3nf}^ + Y,{yakCakdl + h.C.^ (1) 

ak a B ak 
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consisting of three parts H = H^ath + Hioc + Hhyf,. Bath fermionic operators 
are denoted by c^^, cj!^^ while da, da stand for the impurity operators and 
Ha = dtda are the corresponding occupation number operators. The flavor 
index a represents both the spin and the orbital quantum number. Our aim 
is to evaluate the Green's function 

Ga{T',T") = -{Td^{T")dUT')). (2) 

We assume the bath and the Green's function to be diagonal in the flavor 
indices throughout the paper but generalization to non-diagonal baths is 
possible. 

The CT-HYB algorithm is based on the expansion of the partition func- 
tion in the hybridization strength 

Z = Ti [e^^^] = Tr ^^^^^*>ath+Hioc)j'Q- Hhyb{'r)dT _ 

= ^(-1)'= / dn... dr,Tr[e-^(^--+^-=)//,,,(rO . . . i/,,,(ri)]. (3) 

fc=0 JTk-l 

For Hioc diagonal in the occupation number basis, expansion ([3]) can be 
viewed as a sum over so called segment configurations (^-configurations). In 
a schematic way, it can be written as 

Z = Zbath ^ dr Z{t) = Zbath ^ dr Zioc{t) Zhyb{r), (4) 

The partition function of the bath Zbath presents only an irrelevant multi- 
plication factor. The abbreviation r = {{T/,rf}a} is used here to describe 
all the integration variables, that is a set of start-times and end-times for 
all flavors. We stress that we use symbol r to describe an arbitrary set of 
numbers rf and from interval (0,/3). If r does not represent an allowed 
^-configuration, Z(t) is zero. Otherwise, we have 

Zioc = exp - ^ UaE^ - UapOaP , 
\ "</3 / 

Zhyb = U^et Fair,' -tJ). (5) 

a 

where and Oap stand for the total length of segments of flavor a and their 
overlap with segments of flavor (3, and Fa(T) is the hybridization function. 
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A random series of Z-configurations is generated using the Metropolis im- 
portance sampling with the probability density Z{t). This series is then 
used to obtain various quantities of interest such as the Green's function (|2]). 
In the standard approach, from each Z-configuration with K segments of 
flavor a one obtains an estimator of the Green's function that consist of K'^ 
delta-functions 

Ga{r', r") = hMt,6{T' - r/)5(r" - r/)) . (6) 

\ »J ' MC 

Here, M^j stands for the matrix element of the inverse hybridization matrix 
of flavor a. With this approach, however, difficulties are encountered when 
the mean perturbation order of the simulation is low. In such a case, the 
estimator consists of only a few 5-functions and a poor statistics results. 
Moreover, due to the discrete character of estimators, a significant statistical 
noise is observed in Ga{iu}n), and even worse in the selfenergy, at higher 
frequencies. 

Recently, a substantial progress has been achieved in fixing these prob- 
lems. The statistical noise in the selfenergy (and the two-particle Green's 
function) can be significantly suppressed when improved estimators based 
on the equation of motion are used [isj. Besides that, filtering out the 
stochastic noise using orthogonal polynomial representation can be used to 
further improve the results [l9|]. Nevertheless, neither of these methods can 
fully avoid the poor statistics at low perturbation orders. 

3. Measurement of the Green's function 

There is not a unique way to estimate the Green's function from the 
random series of Z-configurations. In fact, there is a substantial freedom in 
this procedure and expression (Q presents only one of the possibilities. In 
this section, we review the general logic of the measurement and establish 
notation used in the rest of the paper. 

Perturbation expansion of the Green's function can be, similarly to series 
03]), written as 

Ga{r', t") = ^fdT G,(r', r"; r) = ^ f dr GUr',r"\ T)G',,,(r). (7) 

It can be represented as a sum over segment configurations with one pair 
of special start-time and end-time at r' and r" (G-configurations). In CT- 
HYB, series ([7]) is integrated by the Monte-Carlo method with Z[t) used as 
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the probability density for the importance samphng. Unless r' = r", there 
are regions of the r-space where Z{t) is identically zero while Ga{T' ,t";t) 
is not. Therefore, Z{t) can not be used to directly sample expansion ([7]). 
In the standard approach, to each G-configuration {t',t",ti} we assign one 
Z-configuration T2 so that if Gaij' ,t"; ti) is non-zero, Z{t2) is non-zero as 
well. Then, once T2 is visited during the random walk, Gair'^r"; ti)/Z(t2) 
is accumulated as the estimator for Gair' ,t"). Drawback of this approach is 
that we obtain information about Ga{T',T") only when r' and t" correspond 
to one start-time and one end-time in the visited Z- configuration. 

The standard measurement can be generalized when to each Gai^', t"; ti) 
we assign a set of configurations Z{t2) with some weight distribution w{t\ r", ti, T2) 
and accumulate 

f, „ Gc,{t',t";ti) „ .Gioc{r',T",Ti)Ghyb{ri) 

^(,T2j ^loc\T2)^hyb\T2) 

(8) 

when any T2 is visited. In principle, the weight function can be chosen 
arbitrarily as long as normalization condition 



J dT2w(r',T",Ti,T2) = l (9) 



is fulfilled for each r', r" and ti. However, for the sake of efficiency, w{t', t", ti, T2) 
should be non-zero only when it is cheap to evaluate ratio Ga^r', r"; ti)/Z(t2) 
in equation ([H]). 



4. Improved measurement 

Here, we introduce a choice of the weight function w(t' ,t" ,ti,T2) which 
constitutes the improved measurement. For a given Z-configuration, we as- 
sign a non-zero weight to a set of G-configurations generated by removing of 
a hybridization line followed by a shift of the lone (ia(r"), (il(rj) operators 
to all positions consistent with the remaining segments. This way, from any 
Z-configuration (apart from the zeroth order of the perturbation theory) we 
obtain a non-zero contribution to Ga{T', t") for all values of arguments so the 
estimator of the Green's function is a continuous function of the imaginary 
time. 

In the following, we distinguish two types of G-configurations. Those 
where r' can be shifted to t" without crossing a segment boundary will 
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<(r';T) d^{T";T) 






Figure 1: Relation between the separated G-configuration (top) and a Z-configuration 
(bottom). When we remove the hybridization hne connecting and rj, we can reach 
G-configurations with r" € {rf^iTrf) and t' e (Tj,Tj_^i). The piecewise-exponential r- 
dependence of functions di{T; t) and (i| ^(r; x) is schematicahy depicted in the two graphs 
in the upper part of the figure. Regions with different atomic states are separated by the 
vertical dashed lines. 



be called connected. Remaining ones will be called separated. The Green's 
function can be split into two terms corresponding to sums of series over 
the connected and the separated configurations as 

G«(r',r") = G^(r',r") + G'f(r',r"). (10) 
Their evaluation is discussed in next two subsections. 

4.1. Separated configurations 

First, we discuss the evaluation of G^{r',T") which is somewhat simpler. 
In expression the ratio of the bath traces equals the element of the 
inverse hybridization matrix as in the standard formula ([6]). The local traces 
differ by the position of one start-time and one end-time. At a fixed segment 
configuration, the local operators of flavor a "feel" an effective r-dependent 
potential 

EtiT;T) = Et+Y,Ua^Mr-^r) (11) 
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consisting of the site energy and the interaction contribution due to all 
segments of other flavors. We define a function 

(i„(r; r) = exp JJ daEi{a; r)) (12) 

so that shifting of the creation operator from r/ to r' results in scaling the lo- 
cal trace by factor da(r^] T)/da(T'] r). When annihilation operator is moved 
from Tj to t", the local trace is scaled by da(T";T)/da{Tj;T). The lower 
bound of the integral in equation (fT2!) can be chosen arbitrarily since it influ- 
ences only an irrelevant multiplicative factor that cancels out. The estimator 
for the separated Green's function then reads 

Ga(i-,7-;T)= 2^ Mtj-r-r-e — vJT^ — r^ir ,t ,t,' ,t^;t). 13 
Ij ' dc,{T^;T)da{r';T) 

The summation does not include neighboring pairs of start- and end-times 
since they contribute to the connected part of the Green's function. We 
switched from abstract notation for weight w{r' ,t" ,ti,T2) and we explicitly 
indicate relevant imaginary-time arguments. In this notation, the normal- 
ization condition reads 

n jy dT^w{T',r'\Tt,T'r-T) = l. (14) 

So far, only the support of w was specified while its explicit functional 
form remains undetermined. In the following, we restrict ourselves to a 
factorized form 

w{t', r", r/, r/; r) = w,{t' , r/; r)t/;,(r", rj; r). (15) 

The standard measurement formula (jS]) is restored with the 5-function weight 
Wc = S{t" - r/), Wa - 5{t' - Tj). Since we want to improve on this and 
use the entire accessible imaginary-time interval, the uniform weight Wc - 
l/(rf - rti + mrL, - rf )), - l/irj^, - r/ + /3^(r; - rj^,)) appears like 
a natural choice. However, our numerical experiments showed that it yields 
strong noise especially for strongly interacting systems and low temperatures 
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and it is always outperformed by the "normalization" weight 



Here, Xia,b)(,T) stands for the characteristic function of interval (a,b). With 
this weight, each (5-function of the standard formula is smeared to a 
piecewise-exponential function with unit integral regardless of the initial po- 
sition of r/ and rj. This property makes the algorithm stable at any regime. 

In principle, we can store estimator (fT3|) in any basis during the simula- 
tion. Nevertheless, the imaginary time representation is inefficient because 
we must either introduce a large discretization error or use a fine imaginary- 
time grid which requires many numerically costly evaluations of expression 
f lT3|) . Despite its clear advantages, the recentl y p roposed orthogonal poly- 
nomial representation for the Green's function 19| is not well suited for our 
purpose either and we use the Matsubara basis. Its main advantage lies in 
the factorization property exp(za;„(r' - r")) = exp(ia;„r') exp(-zaj„r") that 
reduces the computational complexity of the measurement. Besides that, it 
is simpler to Fourier transform a piecewise-exponential function than calcu- 
late its Legendre polynomial representation. In Matsubara frequencies with 
normalization weight (fT6|) we obtain 

Gf(^a;„) = 4 T T 'ir'dr"exp(2a;„(r" - r')) (r', r") 
Jo Jo 

= 1 ^ Mt^Cr(^UJn)A'^i^Un)^ (17) 
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where functions A'^^iUn) and C°'{iUn) are given by 

^A'^-)- ^l. rfr ^^^^^_^^wa{r,,r ,r)exp(^a;„r )- ^^^^^^ , 

A'^ilOOn) = \ r'"dT"da{T";T)exp(zUnT"), 

C^duj ) - - r^dr'^^^^^^w (r^ T'-T)ex^(-zu r') - ^^^'"^"^ 
{tUn) - ^ dr ^^^^^_^^w,{T,,r,T)exp{ luj^r ) - ^^^^^^ , 

CriiuJn) = ^ r dT'd-J{r';T)exp{-iUnT'). (18) 

In principle, the piecewise-exponential functions da{T'; r), da^r"; r) could 
be modulated by a to some extent arbitrary function if a r', T"-dependent 
weight was chosen. However, the normalization condition f[T^ does not allow 
to chose the modulation independently for each t^^ and this makes the 
choice of a better weight non-trivial. It remains an open question whether 
one can choose it in a more clever way so that it will, for instance, simplify 
numerical evaluation of equation f|T8|) . 



4-2. Connected configurations 

There are two different simple ways to reach a connected G-configuration 
from a Z-configuration. The first option is to remove a hybridization line 
spanning a single segment or antisegment and shift da{T"), (io(r') to a new 
position. The second option is to insert a da{T"), d^r') pair into a Z- 
configuration. This approach is routinely used for estimation of the occu- 
pation numbers, that is the equal time Green's function. It can be used for 
measurement of the full Green's function for Hamiltonians including the spin- 
flip terms [oj . But for density-density interactions, the insertion measurement 
alone is not ergodic since it can not access the separated G-configurations. 
The two possibilities are shown in Figure [2] and we will discuss them sepa- 
rately. 

4-2.1. "Remove and shift" approach 

From the two types of G-configurations Gh and Gc depicted in Figure |2l 
we use a non-zero weight only for the later ones. Then, estimator for G'^{iUn) 
is derived analogously to equation ( IT3l) . The only difference here is that the 
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Gc 



Z 
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Figure 2: Different ways how a connected Green's function configuration can be reached 
from partition function configuration Z . Ga corresponds to insertion method. Gb and Gc 
can be obtained by removing of a hybridization line and shifting the remaining operators. 



sum on the right hand side goes over complementary indices and t' and t" 
must be properly ordered. In Matsubara frequencies we obtain 



(19) 



' MC 



where 



^"(^^-) = ^ r^d^' ^\ ^\ ^i^t.rlr\ r"- r) exp(.a;.(r- - /)), 



,(t^;t) da{T';T) 



Again, we can ensure the normalization with a weight function 



w{t:,t^,t',t";t) = N 



(20) 



(21) 



where the multiplication factor is chosen so that condition (Q is fulfilled. 
In the definition of P^{iujn) we have 



and analogously for H^{iUn)- 



(22) 
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Individual estimators given by equation ( IT^ do not have the exact l/o; 
asymptotic. Therefore, the stochastic noise in the Green's function decays 
with the inverse frequency and it is correlated. However, when apart from 
the Green's function we accumulate also the leading order coefficient in 
G^{iUn) ^ -iG'^/un, we can at the end of the simulation evaluate the Green's 
function as G'^{iUn)lG'^ . After this correction, the noise decays with the in- 
verse square of frequency and its correlations are suppressed. 

4-2.2. "Insertion" approach 

With the insertion approach, there is only one Z-configuration from which 
a particular G-configuration can be created. Therefore, as in the standard 
measurement, the weight is uniquely determined by normalization condition 
(Q. When there is at least one segment of flavor a, the estimator reads 

G'^ilOOn) = (E ^"(^^n) + m^^A . (23) 
\ i / MC 

Unlike in equation f lT9|) . elements of the inverse hybridization matrix are miss- 
ing here and definitions of P"{iUn) and H°'(iu)n) cover shorter imaginary- 
time intervals (corresponding to only one segment or antisegment) 

^"(^^") = 4 rdT' r'rfr"Ml^exp(za;„(r"-r')), 

^r(zc.„) = -^ r'^'dr' r^^dr"§4^exp(za;„(r"-r')). (24) 

Unlike with the RS method, we obtain non-zero estimator I^{iUn) even 
from Z-configurations with no segments of fiavor a. Then, there are only 
two possible states: the empty-line and the full-line. Given the configuration 
of other segments, probabilities of the two are given by po = 1/(1 + da{(3]T) 
and pi = da((3; t)/(1 + da((3; r). We can use this and write 

/o(^^n) = ^ . , , / dr' dr" ;) ;^^ xp(za;„(r"-r-)). 

P'^ 1 + da{p]T) Jo Jt' da{T']T) 

(25) 

This estimator is accumulated regardless of whether we reach the empty-line 
or the full-line state. 

With the insertion method, each estimator obeys the exact Ga{iujn) 
l/ioJn asymptotic because the sum of lengths of all segments and antisegments 
equals to /3 and estimator ( !25l) is normalized by the probability prefactor. 
Therefore, the stochastic noise scales with w"^ at high frequencies. 
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4-2.3. Comparison with RS method 

When the mean perturbation order is very low, the insertion method is 
more accurate then the RS method since it works directly from the zeroth 
perturbation order. However, in most practical situations, it has ergodicity 
issues beyond the weak-coupling regime. In strongly interacting systems, 
segments of different flavors are typically anti-correlated so insertion method 
leads to inefficient importance sampling worsening exponentially with U{t" - 
r')- 

The best performance is delivered when we use a combination of the two 
approaches. At zeroth order, we use the insertion measurement since in this 
case ergodicity is assured by the probability factor in (1251) . For higher-order 
G-configurations we use the RS method. With this combination, there is no 
instability and poor statistics can occur only in the unlikely situation with 
the perturbation order sharply peaked around one. Further, when we refer 
to the RS method, this combination is implicitly assumed. 

4.3. Self energy 

Recen tly, an efficient estimator for the selfenergy in CT-HYB was pro- 



posed in [18|]. The method is based on the observation that the selfenergy 



can be expressed from the equation of motion for the Green's function as 

Fo,{iuJn) , . 

Sa(^Wn) = ^ f. X (26) 

where function Fa{iu)n) is the Fourier transform of the correlation function 
Fair'-r") = - E U^,{Td^(T')dUr")n,(T")). (27) 

This method is superior to direct use of the Dyson equation since the selfen- 
ergy is calculated as a ratio of two quantities and only relative errors propa- 
gate. It can be naturally combined with the present measurement procedure. 
All we need to do in order to use estimator fl26|) is to measure the correlation 
function F^ij' - r"). From the definition of the effective time-dependent 
atomic energy ([TT]) we see that 

ZUa^npir;T) = Ei{T;T)-Ei (28) 

Therefore, regardless of what estimator Ga(T' ,t";t) we use for the Green's 
function, the estimator for the F-function is given simply by 

F„(r', r"; r) = G,(t', r"; T)(i?^(r'; r) - E^). (29) 
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In comparison to measurement of the Green's function itself, this step is 
numerically cheap and requires only an acceptable overhead. 



4-3.1. Dynamic part of the self energy 

The selfenergy can be analytically continued to the real axis using the 



maximum-entropy method [20|, |21| in the similar way as the Green's func- 
tion. However, for the MAX-ENT method it is necessary that the continued 
function converges to zero for large frequencies. The full selfenergy does not 
fulfill this condition because of the constant Hartree-Fock term. Nevertheless, 
its dynamic part 

Sf(za;„) = S,(za;„)-Sr (30) 

does. So, in order to use the MAX-ENT analjd^ical continuation, we have 
to measure the Hartree-Fock selfenergy separately and at the end of the 
simulation, we subtract it from the full selfenergy given by equation fl2B]) . 
We can further reduce error bars of the dynamic part when the estimator 
for E^^ is chosen so that it is correlated with asymptotic behavior of the 
full selfenergy. This is achieved by separate accumulating of the asymptotic 
leading order coefficients G"^ and and using their ratio as the estimator 
for 

^"^'""^-G.(zu;.) ^^^^ 

The values of and F^ depend only on the connected part of the Green's 
function and the F function 



p Jo 



F- = ^f^ dTF^ir„r)-F^iT^,r). (32) 

For the insertion method, these equations reduce to simple expressions = 
1 and F^ = E/3*a f^?^('^/3) while for the remove- and- shift approach no such 
simplifications arise. 



5. Implementation notes 

In this section we summarize all steps necessary for implementation of 
the algorithm described above. We can provide full source code of our im- 
plementation upon email request. 
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First, at the beginning of the simulation, in order to avoid repeated com- 
putations of eigenvalues of the atomic Hamiltonian Eioc, we can store them 
into an array indexed by integer representing the local state (individual bits 
of the index can naturally represent occupation numbers of different flavors). 

Next, given a Z configuration from the Metropolis random walk, we must 
divide the imaginary-time axis to regions with fixed atomic states. Each 
region is bounded by two hybridization events at and r^+i and has length 
Arfc = Tfc+i - Tfc. For each and each Matsubara frequency we compute 
exp(iUnTk) and store it (the factorization property of exponential must be 
used in order to carry out this step efficiently). 

Further, for each flavor a we calculate the effective local atomic energy 
E^j^ from equation (ITTi) . exponential factors exp(-£'^^ATfc) and da{Tk;T) - 
nfjQ^exp(-i?^,Ari) defined by equation f|T2|) . 

For each Matsubara frequency and each region (using the precomputed 
arrays), we calculate the following four quantities 

r-Tk+i 

al{iujn) = / (iTexp(-E„^(r-Tfc))exp(2w„r) = 

■JTk 

(exp(-£;^j.Arfc)exp(iw„rfc+i) - exp(iw„rfc))/(za;„ - E^^), 

rrk+i 

cl{iujn) = / drexp(E„j;.(r-rfc))exp(-2a;„r) = 

Jrk 

{exp^-iUnTk) -exp{E^,^ATk)exp(-iUnTk+i))/{iuJn - -E^fc), 
Pki^u^n) = j^""rfr' j^""rfr''exp(i?i(r'-r''))exp(za;„(r''-r')), 
(exp((^a;„ - i?i)ATfc) - 1 - {tu,, - Ei,)An)/{tUr. - E^.f = 

(a'^{iUn) exp(-iUnTk) - ATk)/{iUn - E^i.), 

hki^^n) = -J, dr' fj^ dr" expiEi,ir' - r")) exp(za;„,(r" - r')) = 

(1 - exp(-(^a;„ - Ei,)An) - (tu,, - E^) ATfc)/(zu;„ - E^.f = 
= {c'^{iun)exp{iujnrk) - Ark)/{iuJn- E^^). (33) 

From these quantities we build-up the functions A'^(iUn), C°'(icOn), P°'{iu!n) 
and H'^iiojn) defined by equations f|T8|) and fl20|) . In order to use the weight 
given iDy equations f|T6l) and fl2Tl) . we also need to evaluate the correspond- 
ing normalization integrals that can be expressed in terms of zero-frequency 
limits aj(0+), cj(0+), p°(0+) and /ij(0+). Unless \E'^k^Tk\ « 1, we can cal- 
culate these values simply by putting iuj^ = in equations fl55]) . Otherwise, 
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we must use the Taylor expansions 

1 



vm = ^ - lE^k^ri + . . . 

W = ^ + ^^iAr^... (34) 



Finally, we have 



^iiii^n) = E dairk;T)a'^(iuJn), where : Tk, = r/, Tk,+i = r/+i, 



Afiicon) 



1 A^jlUn) 



1 

C^(i^n) = ^ E d'^in'^'^yti^^n), where : = Tf_i, Tfe^+i = rf , 



Cf(ia;„) 



1 C^jliOn) 

(3 C^O,) ' 



where : r^^ = r^^, r^^+i = r/+i, 

1 



k=ks k=ksl=k+l (^a{Ti,T) 

where : Tk, = r/, r^^+i = rf+i, 

^^'""^ - H^iO,) ■ ^^^^ 

Evaluation of the double sums in definitions of P^{iujn) and Hf{iujn) can 
be accelerated when results from evaluation of A'^{iojn) ^-nd C'^{ioOn) are 
reused and its cost is linear in number of covered regions kg - kg- The zeroth 



15 



order term fl25l) can be evaluated analogously, but order in which individual 
summations are carried out must be chosen based on sign of da{/3;T) - 1 in 
order to avoid the numerical instability. Finally, with results of equations 
fl35|) . we can compute estimators f|T7|) and f|T9|) . 

5.1. Computational complexity 

The obvious disadvantage of the algorithm described in this paper com- 
pared to the standard measurement is the increased computational cost. 
Unlike the standard approach, our measurement is causing a non-negligible 
slow-down of the entire simulation. 

Overall, the computational complexity of the improved measurement 
scales quadratically with the perturbation order and linearly with the num- 
ber of measured Matsubara frequencies because of evaluation of equation 
f|T7|) . Nevertheless, the standard estimator ([6]) has a similar form and it is 
not the real computational bottleneck for any reasonable perturbation or- 
der. In practice, the most expensive part of our algorithm is computation of 
the functions A^{iUn), C^^iujn), P^'i^iuJn) and H°'{iujn) from equation fl5Sl) 
(including all of its inputs). Its cost, however, is only linear in the total 
number of regions, Matsubara frequencies and flavors. Therefore, relative 
cost of measurement compared to one MC step is decreasing with increasing 
perturbation order (temperature). Nevertheless, it is essential not to perform 
measurement after every Monte- Carlo step but only approximately once the 
autocorrelation time to keep the simulation speed reasonable. 



6. Results 

In this section, we present the numerical results for two test cases. First, 
the single band Anderson impurity model with the semielliptic density of 
states (half band- width is set to one). Second, as an example of a realistic 
material calculation, we show results for the cobalt atom in LaCoOs. If not 
explicitly stated otherwise, we use the RS measurement. As a benchmark, 
results obtained with delta estimators measured to Matsubara frequencies 
and Legendre polynomials are used. Our implementation is based on the 



hybridization-expansion solver [22[ from the ALPS library version 2 |23l . |24 
All the other results shown here were also obtained using the ALPS library. 

In figure [3], the impurity selfenergy for the half-filled Anderson model for 
U = 1 and /3 = 10 (mean perturbation order 1.8) is plotted. The results were 
obtained from the equation of motion f l26|) using independent simulations and 
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Figure 3: The real and the imaginary part of the selfenergy for J7 = 1 and /3 = 10. The red 
and the green hne were obtained from delta estimators measured into Matsubara frequen- 
cies and Legendre polynomials. Inset in the right graph shows the difference between the 
measured selfenergy and the analytically computed high-frequency tail. 

the same amount of CPU time. We intentionally use short MC runs here in 
order to make the statistical noise visible to the naked eye. While noise from 
the standard method is increasing with frequency, no such effect is visible 
when continuous estimators or measurement into orthogonal polynomials 
are used. With the orthogonal polynomials, however, the stochastic noise is 
suppressed at the price of introducing strong statistical correlations between 
the values at different Matsubara frequencies. Correlations in the imaginary 
part of the selfenergy obtained from the RS method are much weaker. This is 
demonstrated in the inset which shows the difference between the measured 
selfenergy and the analytically computed high-frequency tail. Errors of the 
real part of the selfenergy are strongly correlated because of uncertainty in 
the Hartree-Fock term. These correlations can be eliminated in the dynamic 
selfenergy by means of equation ( IHT]) . 

Next, we show comparison of accuracies of different methods in figure 
m We plot the relative error, that is, the ratio of the errors obtained with 
continuous estimators and the error from the standard formula. The com- 
parison between the insertion method and the RS method shows that while 
the first one has better high-frequency behavior, it is (exponentially) unsta- 
ble in the strong coupling regime. On the contrary, the RS method yields 
more accurate results regardless of U. The temperature dependence of the 
relative error from the RS measurement shows that benefits of this approach 
persist even in higher perturbation orders. The correction on the exact l/u 
asymptotic discussed in section 14.2.11 further increases the accuracy. 

In figure O the selfenergy of lanthanum cobaltite is plotted. The impurity 
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Figure 4: Left: Comparison of relative error bars of the imaginary part of the Green's 
function for U = 1,2,3,4 and 5 and /3 = 5. Errors from the insertion approach (labeled 
with triangles) are lower in the weak-coupling regime, but they grow exponentially with 
increasing interaction. On the contrary, the RS measurement (labeled with circles) yields 
almost interaction-independent results. Right: Relative errors for [/ = 1 and different tem- 
peratures with (squares) and without (triangles) the correction on the correct asymptotic 
behavior. The mean perturbation order here grows approximately linearly with tempera- 
ture up to {K) « 16 at /3 = 80. 

model used describes the full d-shell of cobalt with five orbitals. The per- 
turbation order is strongly imbalanced here. While the mean perturbation 
order of t2g orbitals is very low, {K) r; 0.3, orbitals are strongly hybridized 
with {K) !« 42. The Cg selfenergy measured with the RS method is less 
accurate than the one obtained from the standard formula apart from the 
high-frequency tail. The reason for this is that the same amount of CPU 
time was used for both simulations and only approximately 60% of Monte 
Carlo steps were done with the slower RS method. 

7. Conclusions and outlook 

The presented measurement algorithm for the Green's function has sev- 
eral possible apphcations. First of all, it solves the poor statistics problem at 
low perturbation orders. At higher orders, it yields stochastic noise decreas- 
ing with Matsubara frequency so it can generate better data for the numerical 
analytical continuation of both the Green's function and the selfenergy. In 
principle, it also allows for measurement of a general local susceptibility of 
non density-density type that is very difficult to obtain with the standard 
method. 

Its main disadvantage is relatively high computational complexity result- 
ing in slowdown of the entire simulation and slightly worse accuracy at low 
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Figure 5: Top: The selfenergy of lanthanum cobaltite calculated with delta estimators and 
the RS method. The two branches correspond to the weakly hybridized t2g orbitals and 
the strongly hybridized eg orbitals. Bottom: Relative errors of the full selfenergy and its 
dynamic part. 
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frequencies. Moreover, since our method essentially relies on fast and simple 
evaluation of the local trace, generalization to the matrix formalism would 
be difficult, if possible at all. 
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